Computational framework for the generation of one-dimensional vascular models accounting for uncertainty in networks extracted from medical images

One-dimensional (1D) cardiovascular models offer a non-invasive method to answer medical questions, including predictions of wave-reflection, shear stress, functional flow reserve, vascular resistance, and compliance. This model type can predict patient-specific outcomes by solving 1D fluid dynamics equations in geometric networks extracted from medical images. However, the inherent uncertainty in in-vivo imaging introduces variability in network size and vessel dimensions, affecting hemodynamic predictions. Understanding the influence of variation in image-derived properties is essential to assess the fidelity of model predictions. Numerous programs exist to render three-dimensional surfaces and construct vessel centerlines. Still, there is no exact way to generate vascular trees from the centerlines while accounting for uncertainty in data. This study introduces an innovative framework employing statistical change point analysis to generate labeled trees that encode vessel dimensions and their associated uncertainty from medical images. To test this framework, we explore the impact of uncertainty in 1D hemodynamic predictions in a systemic and pulmonary arterial network. Simulations explore hemodynamic variations resulting from changes in vessel dimensions and segmentation; the latter is achieved by analyzing multiple segmentations of the same images. Results demonstrate the importance of accurately defining vessel radii and lengths when generating high-fidelity patient-specific hemodynamics models.


INTRODUCTION
Medical images, including computed tomography (CT) and magne:c resonance imaging (MRI), are widely used for disease diagnos:cs and treatment planning (Chassidim et al., 2015;Doi, 2006;Khoo et al., 1997;Alas: et al., 2006).They are non-invasive and reliable ways to assess the structure of the heart, lungs, and vasculature.However, such images do not provide insight into hemodynamics, which is essen:al to determine the health of the cardiovascular system.One way to assess cardiovascular health is to embed onedimensional (1D) pa:ent-specific computa:onal fluid dynamics (CFD) models within the image.This can be done by extrac:ng a pa:ent-specific vascular network from the medical image, and then solving fluid dynamics equa:ons in each vessel within the network to make predic:ons of dynamic blood pressure and flow.To generate high-fidelity predic:ons for clinical applica:ons, it is essen:al to minimize and quan:fy the uncertainty associated with this process.In this context, uncertainty arises from two main sources: the network extracted from the medical image and parameters in the computa:onal model.Numerous studies have explored how varia:ons in model parameters influence fluid dynamics predic:ons (Bertaglia et al., 2020;Chen et al., 2013;Ninos et al., 2021;Ye et al., 2022), but only a few (Colebank et al., 2019;Sankaran et al., 2015) have addressed the role of varia:ons in network geometry.This study focuses on determining varia:ons in vessel radius, length, and network connec:vity and exploring how it affects 1D hemodynamic predic:ons in arterial networks.
Labeled trees are generated from three-dimensional (3D) rendered surfaces obtained using image segmenta:on.Image segmenta:on allows for analysis of anatomical data by isola:ng key features in an image, in this case vascular networks (Pham et al., 2000;Pa:l and Deore, 2013).However, there is inherent uncertainty in the image segmenta:on process due to complex anatomy, pa:ent mo:on during scanning, variability in user exper:se, and noise (Sharma and Aggarwal, 2010;O'Donnell, 2001;van Rikxoort and van Ginneken, 2013;Buelow et al., 2005).Despite advances including semi-and fully automa:c extrac:on using machine learning, most soUware requires manual edi:ng to capture the en:rety of the desired anatomical features (van Rikxoort and van Ginneken, 2013).In par:cular, pa:ent images are only captured up to a finite resolu:on, and certain anatomies, such as the rapidly branching structure of the pulmonary vasculature, are nearly impossible to extract in an automated manner (Colebank et al., 2019;van Rikxoort and van Ginneken, 2013).Other anatomies, such as the aorta, are difficult to automa:cally extract due to differences in spa:al scales within the structures we wish to capture, which includes the smaller branching neck vessels and larger main aor:c branch.Overall, a major challenge is that vascular networks span several orders of magnitude, but resolu:on within the image is constant (Lakhani, 2020), adding significant uncertainty when extrac:ng data, especially for small vessels (Colebank et al., 2019).Moreover, for images acquired without contrast, differen:a:ng arteries, veins, and other structures is not trivial (van Rikxoort and van Ginneken, 2013).This study generates networks from images segmented using 3D Slicer developed by Kitware, Inc. (Federov et al., 2012), which is a widely used open-source soUware.
AUer segmen:ng the image and genera:ng a 3D rendered surface, centerlines are iden:fied along each vessel, crea:ng a 1D characteriza:on of blood vessel geometry (Pham et al., 2000;Pa:l and Deore, 2013).This study builds networks using the Vascular Modeling Toolkit (VMTK) (Izzo et al., 2018), which generates centerlines by placing maximally inscribed spheres along each vessel.This process is sensi:ve to intricate vascular networks containing vessels with high tortuosity, rapid branching structures, and narrow lumen.In these networks, maximally inscribed spheres placed by VMTK intersect earlier than expected, leading to junc:on nodes placed outside of physiological domain.In par:cular, VMTK oUen posi:ons junc:on nodes prior to the os#um region, the area surrounding the opening of a blood vessel where it divides into mul:ple vessels and where the radius is not defined (Cheng, 2019).This limita:on has been described in several studies, including Ellwein et al. (2016) and Pfaller et al. (2022).To address the limita:on, Ellwein et al. (2016) adjust junc:ons by examining the angles between vessels, while Pfaller et al. (2022) u:lize slices rather than maximally inscribed spheres to generate centerlines.Other programs exist to generate centerlines, such as SimVascular (Updegrove et al., 2017), CRIMSON (CardiovasculaR Integrated Modelling and Simula:ON) (Arthurs et al., 2021), and SGEXT (Cappek et al., 2022), each balancing accuracy and efficiency.Independent of the soUware used to generate centerlines, the vessel radius is not defined in the os:um region, yet fluid dynamics models solving the equa:ons in 1D trees require a radius measurement at all points along the vessel.To address this limita:on, this study devises an algorithm to move junc:ons into the os:um region and uses sta:s:cal change points to iden:fy a segment within each vessel's centerlines that best represents its radius.
CFD modeling provides insights into the func:on of the vasculature by predic:ng hemodynamic proper:es that cannot be measured in vivo (MiMal et al., 2016;Chinnaiyan et al., 2017).By integra:ng with pa:entspecific networks, these methods can be used to study the outcome of treatments or predict disease progression (Candreva et al., 2022).Several recent studies have used pa:ent-specific geometry derived from medical images to predict hemodynamics (Taylor-LaPole et al., 2023;Colebank et al., 2019;Bartolo et al., 2022;Baksta et al., 2016;Formaggia et al., 2006;Gray and Pathmanathan, 2018;Yang et al., 2019).They superimpose dynamic pressure and or flow data onto geometric domains extracted from images and make predic:ons of blood pressure, flow, shear stress, and wave intensity to simulate the impact of disease and treatment.The most detailed models include 3D CFD and fluid-structure interac:on, which provide predic:ons of local flow paMerns important in regions with significant secondary flows.Unfortunately, 3D models are computa:onally expensive and challenging to calibrate to pa:ent data (Yang et al., 2019;Botnar et al., 2000).Thus, this study focuses on 1D fluid dynamics models, providing an efficient alterna:ve to predict volumetric flow and pressure averaged over vessel cross-sec:ons.This model type is useful for predic:ng wave-propaga:on and wave-intensity, organ perfusion, func:onal flow reserve, vascular resistance, and compliance (Shi et al., 2011).In addi:on, 1D models can be easily calibrated to data, making them ideal for pa:ent-specific simula:ons, and provide reasonable agreement to 3D models (Moore et al., 2005;Reymond et al., 2012;Blanco et al., 2018).
Uncertainty arises in nearly every level of the modeling process, including in the assump:ons to derive the equa:ons that model the physical system, parameter choices, and geometry extracted from image segmenta:on.Uncertainty quan:fica:on has been studied extensively in computa:onal modeling of hemodynamics.For example, Colebank et al. (2019) study the influence of pre-segmenta:on parameters, such as thresholding and smoothing, on 1D predic:ons in ex vivo murine arteries.Chen et al. (2013) incorporates intrinsic parametric uncertain:es to depict blood flow and pressure within a stochas:c model of arteries.Bertaglia et al. (2020) highlight the effects of changes in elas:c and viscoelas:c parameters in fluid-structure interac:on models of arteries.While these studies effec:vely quan:fy uncertainty in the modeling and parameter es:ma:on process, they do not examine uncertain:es associated with the network extrac:on process.It is well known that there is a significant interplay between vascular geometry and hemodynamics, and therefore, the uncertainty of geometric features of vessels must be explored and quan:fied (Gounley et al., 2017).
This study inves:gates the process of extrac:ng networks from medical images and the associated hemodynamic uncertainty.We use 3D Slicer (Federov et al., 2012) to generate 3D rendered surfaces from CT images of the pulmonary arteries for a normotensive human subject and from magne:c resonance angiography (MRA) from a double outlet right ventricle (DORV) pa:ent.VMTK is used to generate centerlines along each vessel.From these, we generate labeled directed trees characterized by vessel connec:vity, radius, and length.To improve fidelity of the generated models, we devise algorithms to align networks and adjust junc:ons to physiologically accurate loca:ons.We then employ sta:s:cal change point analysis to determine a representa:ve radius for each vessel, avoiding outlier data.We propagate uncertainty in segmenta:on, junc:on loca:on, and radius to a 1D fluid dynamics model, finding that varia:on in vessel radius have the most significant impact on fluid dynamics predic:ons.In addi:on, we generate mul:ple segmenta:ons from each image -five for the pulmonary vasculature and ten for the aorta -to study varia:ons between segmenta:ons.Finally, in the pulmonary vasculature, each segmenta:on captures a varying number of vessels, so to enable comparison, we standardize networks using a radius pruning algorithm.Our framework generates high-fidelity labeled trees and provides hemodynamic predic:ons with uncertainty within expected measurement error.As pa:ent-specific modeling becomes increasingly popular in medical applica:ons, we argue that quan:fying geometric uncertainty using our methodology (illustrated in Figure 1) is cri:cal to obtaining reliable hemodynamic predic:ons.

Image Analysis
We analyze two images: a chest CT image from a healthy, 67-year-old female volunteer captured using the Omnipaque 350 contrast agent available at the Vascular Model Repository (Wilson et al., 2013) (hMp://simvascular.github.io/)and an MRA image from a DORV pa:ent acquired at the Texas Children's Hospital Heart Center.Data collec:on was approved by the Baylor College of Medicine, Ins:tu:onal Review Board (H-46224: "Four-Dimensional Flow Cardiovascular Magne:c Resonance for the Assessment of Aor:c Arch Proper:es in Single Ventricle Pa:ents").The MRA image was obtained using :me-resolved, contrast-enhanced imaging performed using localizing sequences from the cardiac MRI with a 0.1 mL/kg intravenous gadolinium contrast bolus injec:on.More details on the image analysis protocol can be found in (Taylor-LaPole et al., 2023).

3D
Rendering.The pulmonary and aor:c networks are rendered as a 3D surface using the open-source image segmenta:on soUware 3D Slicer developed by Kitware, Inc. (see hMp://www.slicer.org)(Federov et al., 2012;Kikinis et al., 2014).Image intensi:es in the CT image range from 50 − 3027 Hounsfield unit (HU) and from 75 − 327 HU for the MRA.Image voxel size for the CT image is 0.68 × 0.68 × 5.00 mm 3 and is 1.22 × 1.22 × 1.40 mm 3 for the MRA.The CT image was segmented five :mes and the MRA image 10 :mes by the same user, each using the same threshold and smoothing factors.Example 3D rendered surfaces for the two geometries are shown in Figure 2.
The MRA segmenta:on was acquired u:lizing 3D Slicer's thresholding, growing from seeds, cukng, and smoothing tools in a semi-automa:c manner.To segment pulmonary arterial networks from the CT image, we first iden:fy the main (MPA), right (RPA), and leU (LPA) pulmonary arteries through thresholding followed by manual pain:ng, erasing, and cukng.This is followed by manual segmenta:on of the lobar, segmental, and subsegmental vessels (Colebank et al., 2021b).The rendered surfaces are saved as standard tessella:on language (STL) files, which can be opened in Paraview, Kitware Inc. (Ayachit, 2015).The skeletonized images are converted to a VTK polygonal data file and then processed through VMTK (An:ga et al., 2008) (see hMp://www.vmtk.org).
Centerlines.VMTK (Izzo et al., 2018) determines centerlines in the 3D rendered surfaces of arteries using maximally inscribed spheres.Spheres, which touch at least four points on the 3D rendering, are inscribed throughout each vessel, and the medial axis is approximated by a Voronoi diagram (An:ga et al., 2008).Centerlines are defined as the minimal path along the inverse of the spheres' radii (An:ga et al., 2008).User-specified inlet and outlet points establish centerline boundaries, guiding VMTK in a recursive process star:ng at the terminal vessels.When two centerlines intersect, a junc:on node is placed.Example 3D rendered surfaces with centerlines are shown in Figure 2.

Labeled Tree Genera3on
A tree is generated to form arterial networks and used in fluid dynamics simula:ons.This tree consists of two parts: (a) large arteries iden:fied from centerlines derived from MRI and CT images and (b) small vessels represented by fractal trees.The small vessels extend the networks from (a) represen:ng vessels that are not visible in the images.The caliber of vessels represented by fractal trees will depend on the image resolu:on.
Custom MATLAB algorithms are used to construct a raw labeled directed tree from VMTK centerlines (Colebank et al., 2019(Colebank et al., , 2021a)).This algorithm transforms centerlines to edges with $x,y, \text{ and }z$ coordinates at the center of each maximally inscribed sphere and defines junc:ons as nodes where two or more centerlines intersect.The tree has three types of vessels: the root vessel, which starts at the inlet and ends at the first junc:on node, central vessels, which begin and end at junc:on nodes, and terminal vessels, which start at a junc:on node and end at an outlet node.Each vessel contains  nodes with associated , , and  values, radius, and length.The radii along the vessel are determined from the maximally inscribed spheres.Each vessel is assigned a length specified by calcula:ng the Euclidean distances between , , and  coordinates of all nodes within the vessel.Informa:on about each vessel is stored in a connec-:vity matrix that links parent and daughter vessels to describe how vessels are connected in the tree.
The ten aor:c networks all have 13 bifurca:ng vessels joined using the same connec:vity matrix, but vessel length and radii vary between networks.In contrast, the rapidly branching pulmonary network requires manual segmenta:on to separate arteries from veins and occasionally leads to trifurca:ons.Therefore, in addi:on to vessel radius and length varia:on, the five segmented pulmonary arterial networks have different numbers of vessels and connec:vity.The trees range from 167 to 238 vessels with 0.5 − 2% of parents having three daughters and the rest of the parents having two daughters.There are no quadfurca:ons or higher in the networks we segment.To compare these networks, similar to our previous study (Miller et al., 2023) we employ a pruning procedure to generate trees with the same number of vessels.This involves itera:vely removing the smallest terminal vessels with terminal siblings un:l the networks reach the desired number of vessels (Figure 3).In this case, we prune networks un:l they reach the size of the smallest network of 167 vessels.As seen in Figure 3, dis:nct networks have a more similar radius distribu:on aUer pruning.Small Vessels.Due to limits in image resolu:on, it is not possible to obtain measurements for small arteries and arterioles.To represent all arteries in the network, star:ng at the heart and ending at the capillary level (Boron and Boulpaep, 2017), we augment the large vessel domain with asymmetrically-branching, self-similar fractal trees (Olufsen et al., 2000).These structured trees have a root vessel with the same radius as the large terminal vessel it is aMached to.Each parent vessel, , branches into two daughter vessels  = ( !,  " ), with length and radii defined as where ,  < 1,  +  ≥ 1, and ℓ %% is a length to radius ra:o.Values for these parameters are listed in Table 1 and are assigned based on literature and analysis of murine micro-CT images (Olufsen et al., 2000;Chambers et al., 2020).The structured trees bifurcate un:l a specified minimum radius,  &'( , is reached.For this study,  &'( = 0.001 cm, consistent with a red blood cell's diameter (Townsley, 2012).

Adjusted Labeled Trees
The directed trees represen:ng the large vessels have two flaws: (1) the junc:on nodes are oUen placed outside of the os:um, and (2) the low-dimensional representa:on of the vessel radius is incorrectly es:mated if all nodes within a vessel are used.In par:cular, the radii assigned for vessel nodes within the os:um do not represent the vessel's true radius, as it contains values for mul:ple vessels as they split.To mi:gate these limita:ons, we have designed an algorithm that moves the junc:on node to the center of the os:um and iden:fies the segment within each vessel that best represents its radius.Nodes along this op:mal segment are then used to determine a radius and its standard devia:on along the en:re vessel.The input to our algorithm is a connec:vity matrix, vessel length and radii generated by VMTK, and postprocessed by custom algorithms devised in this study.Results include an adjusted labeled tree that reduce the uncertainty in fluid dynamics predic:ons propagated from uncertainty in vessel radii and length.

Adjusted Junc9on Nodes
Step 1: Star:ng at the VMTK-derived junc:on nodes  ) = C ) ,  ) ,  ) D that connect terminal vessels, we compute the distance between the , ,  coordinates of the daughter vessels as where  ,* = ( ,* ,  ,* ,  ,* ) denotes the th point aUer the junc:on node  ) for each of the daughter vessels  = ( !,  " , … ,  # ), where # is the number of daughter vessels and  ,.denotes the last node of the shortest daughter vessel.If there is a trifuca:on, we consider the pair of daughters that has the smallest distance between them.For each node  ,* , we then normalize  * by dividing by the maximum distance between the daughter vessels,  &/0 .A new junc:on node  K ) is placed at the smallest  for which  * / &/0 > 0.1.The cutoff value 0.1 is an empirically chosen input parameter.Figure 4(a) shows the th junc:on node represen:ng the end of the parent vessel that separates into a bifurca:on and nodes along each of the daughter vessels  ,* , where  ! is red and  " is blue. Step Figure 4(b) shows the new parent vessel nodes (green) and the adjusted junc:on node  S ) (purple).Note that the path from the new junc:on to the daughter vessels is not smooth.Step 3: The new junc:on node  S ) is connected to nodes  ,* ,  = 1, … ,  along each daughter vessel.To ensure a smooth transi:on, we compute a weighted average in sec:ons where daughter vessel centerlines are within 10 − 20%, i.e., nodes for which 0.1 <  * / &/0 ≤ 0.2,  = , … , .Here,  is the index of the first node for which  *  &/0 ⁄ ≥ 0.2.The new daughter nodes,  S ,* ,  = 1, … , , shown in Figure 4(c), are computed as a weighted average between the extension of the parent branch (dashed green line on Figure 4(c)) and the original nodes  S ,* ,  = , … , .The extension of the parent branch is computed by averaging all nodes  = , … ,  along the daughter vessels as described in Equa:on (3).
The degree of weigh:ng between the extended parent vessel and the daughter vessels is computed based on the ra:o of the radius and the distance between the extended parent and the daughter vessels.The loca:ons of the new daughter nodes are calculated as Each daughter is individually averaged with the extended parent using weights, i.e,  !,  " , … ,  # are averaged separately.For the pulmonary arteries, we use a double-weighted average ( = 2).Since the aorta consists of a large, curved main vessel, centerlines lie closer to the main aor:c vessel than the neck vessels, branching upward at a nearly 90 −degree angle.For these vessels, we calculate a 100 −:mes weighted average  = 100 towards the main aor:c branch.The weigh:ng cutoff of 0.2 and the degree of weigh:ng  are empirically chosen input parameters that can be adjusted to fit the par:cular applica:on.Based on these adjustments, new lengths and radii are assigned to the labeled directed tree.

Adjusted Radii
The radii iden:fied by VMTK at each node along vessels include values in the os:um region, oscilla:ons associated with irregulari:es on the surface of the 3D rendering, and regions close to image resolu:on.Therefore, to determine the radius of each vessel, we devise a three-step process that (1) uses change points to characterize the vessel shape, (2) determines the region that represents the vessel radius (Algorithms 1 and 2), and (3) checks the consistency of the generated labeled tree. Step

Number of change points.
We use the R/changepoint package (Killick and Eckley, 2014) to determine the op:mal number of change points in each vessel.This program employs binary segmenta:on, itera:vely par::oning the dataset and fikng a piecewise constant line in each segment.To determine if a change has occurred, the algorithm seeks the loca:on at which the maximum of the log-likelihood is aMained, where  g ! and  g " are es:mates of mean and variance before and aUer the data is par::oned.We use Bayesian Informa:on Criterion () to penalize the where  is the number of parameters and  is the number of data points.Here,  = 2 since parameters are the mean and variance.Once the  is minimized, a change point is placed at max We limit the number of change points to three to avoid overfikng and enable consistent analysis between vessels.This binary change point algorithm efficiently iden:fies the number of change points at a low computa:onal cost but at the expense of precise change point loca:on (Eckley et al., 2011;Killick et al., 2012;Killick and Eckley, 2014), so we employ R/segmented to improve change point placement.

Change point loca6on.
Op:mal change point placement is determined using the customizable and flexible R/segmented package, which requires specifica:on of the number of change points (selected as described above).Like R/changepoint, this program employs binary segmenta:on, itera:vely par::oning the dataset.Within each segment, the algorithm maximizes the likelihood between a linear regression model and the data (Muggeo, 2003(Muggeo, , 2008) (Figure 5), using a sta:s:cal model that accounts for the slope, mean, and variance, which are parameters of the regression model.
A two-step process is employed using a regression model between the mean response, [], and parameters, , to determine the loca:on of change points, and then a simpler regression model (using the same response and a predictor) to fit piecewise linear func:ons between change points.The regression model used to determine change point loca:on is defined as where (⋅) is the link func:on for [], ℎC g ;  * D is a parameteriza:on for ,  * ,  = 1, 2, 3 are the change points, () is the predictor,  is the length along the vessel.A first-order Taylor expansion around an ini:al The response and parameters form a piecewise rela:onship with the terms Note that ( −  * ) 6 = ( −  * ) × C( >  * )D where (⋅) is the indicator func:on.(⋅) = 1 when the statement is true, and (⋅) = 0, otherwise.Using this parameteriza:on •  is the slope of the line segment prior to the change point.
•  is the difference in slopes of the lines before and after the change point.
• ( + ) is the slope of the line segment following the change point.
The change point is fixed at  (>) = C −  * (>) D 6 and  (>) = −() >  * (>) and iterated  :mes un:l there is no improvement in the  es:mate.Using radii,  (cm), as the response, and  ,* = (1, … , ) the loca-:on along the vessel, the linear regression model is given by where  < ,  ! are regression coefficients.Step 2: Vessel segment selec9on to adjust vessel radius and determine its uncertainty.This step determines the segment within each vessel that best represents the vessel radius and from which the value can be reliably es:mated.We refer to this as the op:mal radius segment,  ?@A .Given the varia:on in vessel morphology, a mul:step process is used to iden:fy  ?@A .The process is described in detail in Algorithms 1 and 2 with R and MATLAB code available in the repository found at hMps://github.com/msolufse/CDG_NCSU/tree/master/VascularTreeFromImages.
Most 1D fluid dynamics models predict hemodynamics, assuming that vessels are either straight (Colebank et al., 2021a) or tapering (Olufsen et al., 2000).Typically, in rapidly branching vascular networks such as pulmonary vasculature, the brain, or the liver, vessels do not taper significantly between junc:ons, whereas in networks with longer vessels, e.g., networks including the aorta, caro:d, brachial, or iliac arteries, vessels taper along their length (Taylor-LaPole et al., 2023;Olufsen et al., 2000;Epstein et al., 2015).
The fluid dynamics code used in this study allows vessels to be straight or taper exponen:ally.Therefore, the proposed algorithm generates radii labels for both cases.
Non-tapering vessels are assigned a constant radius along the vessel, averaging values in  ?@A , such that where  ?@A () refers to the th radius observa:on in  ?@A and  is the number of points in  ?@A .Note that although  " is assigned a single value based on the mean, it is a vector because we assign this value to every node along the length of the vessel.In tapering vessels, we fit a decaying exponen:al func:on to radius values in  ?@A given by  "() =  '( exp( log( '(  ?CA ⁄ )), where  '( and  ?CA refer to inlet and outlet radius of the vessel.The parameters  '( and  ?CA are es:mated using MATLAB's op:miza:on func:on fminsearch() (MathWorks, 2022) minimizing the least squared error For all vessels, the standard devia:on of the radius is determined by The following sec:ons describe how  ?@A is selected for the two vessel types.
Non-tapering vessels.To iden:fy  ?@A for non-tapering vessels, the slopes of regression lines fit to  * are considered, choosing the sec:on with the smallest absolute value slope,  &'( .The smallest slope is found as  &'( = min( * ),  = 2, … ,  + 1.Note that the first segment (the sec:on prior to the first change point) is not included as it is typically within the os:um region.To account for the varia:on in vessel morphometry, this process requires user specifica:on of three parameters: the minimum percentage needed to ensure the length of single segment has sufficient data (here set at /4 = 25%), the threshold for analysis of mul:ple segments to specify  ?@A (here set to /2 = 50%), and a threshold  to compare the slopes of the last two segments (here set to  = 2.8).
Case 1: Iden:fy the segment with the smallest absolute value slope located between change points  * ,  = 1, … ,  and , the last node in the vessel.For a vessel with three change points, there are three sec:ons to consider, ( !,  " ), ( " ,  E ), or ( E , ).If the sec:on with the smallest absolute value slope includes more than 25% of the total vessel observa:ons, we assign it to be  ?@A .The threshold /4 = 25% is one of the three input parameters.
Case 2: For vessels not included in case 1, we determine the first ( = 1) and last ( = ) radius values in  ?@A () by analyzing change point loca:ons.
(a) To identify the first point in  ?@A , we consider the location of the first change point,  ! .If  ! is position beyond the 25% mark in the vessel,  ?@A begins at ( ! ).If  ! is located before the 25% mark, the first change point is assumed to be in the ostium region, so we start  ?@A at the second change point, ( " ).To avoid excluding data, when  " is located more than 50% into the vessel, we instead start  ?@A at the 25% mark in the vessel, (⌈ 4 ⁄ ⌉).
(b) The last radius value in  ?@A is determined by considering the last change point,  7 , or the last radius observation, ().Which option is used depends on the data structure.We compare the slopes of the second-to-last,  7 , and last,  76! segments as these slopes differ significantly for many vessels, particularly in small vessels with diameters close to the image resolution (Cui et al., 2019;Lesage et al., 2009).To compare these slopes, we compute If  F ≤ , we let  ?@A () = ( 7 ) and otherwise,  ?@A () = ().
Tapering vessels.The methodology for determining  ?@A for tapering vessels is similar.Instead of examining the segment with the minimum slope, we examine the longest segment  &/0 = maxClength( !,  " , … ,  &6! )D .Iden:fica:on of the op:mal radius segment for tapering vessels depends on four parameters: the threshold for fikng  ?@A from a single segment (here 0.4  = 40% of the total vessel length), the parameters used in case 2 for non-tapering vessels (/4 = 25% and /2 = 50%), and the threshold parameter  (for this study,  = 2.8).
Case 1: If the longest segment includes more than 40% of the total vessel length, we let  ?@A =  &/0 .The threshold 0.4  = 40% is one of the four input parameters.
Case 2: For vessels not included in case 1, we determine the first and last points in  ?@A as described in case 2 for non-tapered vessels.
Case 3: For vessels with no change points (typically short vessels), we let  ?@A = (1, ).Note these vessels are assigned a non-tapering radius.
Case 4: Vessels with op:mal radius segment not iden:fied in cases 1-3 , an exponen:al is fit of the same form described above (Equa:on ( 12)).The op:mal inlet radius is determined from (1), and the op:mal outlet radius is determined from ().
Step 2: Network Consistency Check.For healthy subjects, vascular networks sa:sfy that the radius of the parent vessel is larger than the radii of the daughter vessels  $ >   .In most cases, a parent splits into two daughter vessels and very rarely (less than 2%) splits into three (Chambers, 2022).In addi:on, the combined area of the daughters is greater than the area of the parent vessel, such that ∑  >  $ (Zamir, 1978;Murray, 1926;Uylings, 1977;Townsley, 2012).In vessels where these condi:ons are not sa:sfied, we adjust the radii depending on the vessel's posi:on within the network.If one terminal daughter,  !, violates the assump:on, we let  # !=   $ and if two terminal daughters violate the assump:on, we let  # !=   $ and  # " =   $ .Internal vessels or terminal junc:ons with trifurca:ons viola:ng these condi-:ons are manually assigned new radius values from inspec:on of the data.
Algorithm 1: Given radius data ,  ≤ 3 change points , segments between change points , and slopes of lines fitted within each segment , this function identifies the region that best determines the vessel radius,  ?@A , in straight and tapering vessels.

Fluid Dynamics Model
We predict hemodynamics in each vessel using a 1D fluid dynamics model that relates blood pressure, flow, and cross-sec:onal area.The model is separated into two domains predic:ng hemodynamics in (1) large arteries with geometry determined from images and (2) small arteries represented by structured trees.Similar to previous studies (Taylor-LaPole et al., 2023; Bartolo et al., 2022;Olufsen et al., 2000;Colebank et al., 2019), in the large vessels, we solve a non-linear unsteady 1D approxima:on of the Navier-Stokes equa:ons, enforcing conserva:on of mass and momentum, and in the small vessels, we solve a linearized wave equa:on model.Parameters for each domain are listed in Table 1 and the supplement.Large Vessels.Blood pressure (, ) (mmHg), flow (, ) (cm 3 /s), and area deforma:on (, ) (cm 2 ) are predicted under the assump:ons that each vessel can be represented by a deformable axisymmetric cylinder with an impermeable wall.We assume that blood is incompressible, Newtonian, viscous, and homogeneous and that the flow is irrota:onal and laminar.Under these assump:ons, the conserva:on of mass and balance of momentum are given by where µL is viscosity,  =  H / is the kinema:c viscosity,  (cm) is the radius,  (s) denotes the temporal coordinate, and x (cm) is the axial posi:on.We assume a flat velocity profile with a linearly decreasing boundary layer (Nichols et al., 2011;Pedersen et al., 1993;Pries et al., 1992;Olufsen et al., 2000) with thickness  = ¡/2 (cm), where  (s) is the length of the cardiac cycle, i.e., To close the system of equa:ons, we assume that the vessel wall can be modeled as linearly elas:c and isotropic (Qureshi et al., 2019) under which pressure and area can be related as where  (g/cm/s 2 ) is Young's modulus, h (cm) is the vessel wall thickness, p0 (mmHg) is a reference pressure, A0 (cm 2 ) is the cross-sec:onal area, and R0 (cm) the radius of the vessel determined by our algorithm.
At the inlet of the root vessel, we prescribe an inflow waveform extracted from data.At the MPA, we specify a flow waveform using data from SimVascular (Updegrove et al., 2017).In the aor:c vasculature, a flow waveform at the root of the ascending aorta is extracted from 4D-MRI data provided by Baylor College of Medicine and Texas Children's Hospital.At junc:ons, we enforce con:nuity of pressure and conserva:on of flow between the parent vessel p and the k daughters, d = (d1,d2,...,d#)) Small Vessels.In small vessels with radii smaller than image resolu:on, viscous forces are dominant, allowing us to neglect the nonlinear iner:al terms and linearize Equa:ons (15), as done in previous studies (Olufsen et al., 2000).In the frequency domain, the conserva:on of mass and momentum equa:ons become where  < ,  ! are the zeroth and first order Bessel func:ons,  < " =  E  < /,  is vessel compliance, and ℎ/ < has the same form as Equa:on ( 16) but with parameters  !L ,  " L , and  E L .Viscosity becomes more significant for smaller vessels (Pries et al., 1992), so in this domain, we compute viscosity as a func:on of vessel radius as where  = (2 < /(2 < − 1.1)) " is the rela:ve viscosity at a hematocrit level of 0.45.Solu:ons of these equa:ons allow us to compute (0, ) = (0, )/(0, ) at the beginning of each vessel as a func:on of the impedance at the end of the vessel as (0, ) =  P M! sin(/) + (, ) cos(/) cos(/) + _(, )sin (/) .
Using same junc:on condi:ons as the large vessels, the impedance at the root of the structured tree is computed recursively.The combined system of equa:ons is solved numerically using the two-step Lax-Wendroff finite difference scheme (Colebank et al., 2019(Colebank et al., , 2021a;;Bartolo et al., 2022;Taylor-LaPole et al., 2023;Olufsen et al., 2000;Chambers et al., 2020).

Sta3s3cal Analysis and Uncertainty Quan3fica3on
To determine uncertainty associated with image segmenta:on, we perform a one-way analysis of variance (ANOVA), comparing five segmenta:ons of the pulmonary vasculature and ten segmenta:ons of aor:c vasculature prior to adjus:ng the labeled trees (i.e, before modifying junc:ons and extrac:ng radii based on change points).This sta:s:cal test determines if significant differences exist among the variance of different samples (Kauffman and Schering, 2007).We conduct the ANOVA tests in a vessel-specific manner to determine if there are sta:s:cally significant differences in vessel radii and length from the different 3D renderings and centerlines of the same subject.The ANOVA calcula:ons are performed in R using the func:on anova() from the R/stats package (R Core Team, 2013) (Table 2).More details about ANOVA can be found in the supplemental material (S4).
To examine variability within the radius data, we calculate coefficients of varia:on as  =  %̂/ ̂ for vessel radii in each segmenta:on (GraUon et al., 2022).Here, ̂ is the mean of the op:mal radius region and  %̂ is its standard devia:on.
We also compare kernel density es:ma:ons (KDE) of vessel radii for each segmenta:on.The kernel density represents the data using the radii's probability density func:on (PDF).These KDEs are ploMed to visualize the es:mated distribu:on of the radii throughout the segmented vasculatures (Weglarczyk, 2018).We use MATLAB's ksdensity() to convert the radii from each segmenta:on into a PDF (MathWorks, 2022).

Simula3ons
To assess uncertainty in image segmenta:on and 1D network genera:on, we perform two fluid dynamic simula:on types.First, we compare mul:ple networks generated from the same image, assessing the impact of segmenta:on uncertainty and pruning.Then we study the impact of varia:on in junc:on placement and vessel radius within a representa:ve network.
3D Rendering Varia9on.The pulmonary network was segmented five :mes and the aor:c network was segmented ten :mes by the same user.For both, we used the same threshold parameters and smoothing type.All aor:c networks have the same number of vessels and connec:vity, but these features differ in the pulmonary trees.Thus, we examine the effect of discrepancy in vessel count and the impact of pruning.
Varia9ons in Junc9on and Radius Es9ma9ons.A representa:ve network is selected to study the influence of junc:on placement and vessel radius on hemodynamics.We demonstrate the impact of altering the loca:on of the junc:on node, which changes the length and radii considered for each vessel in the system.Once junc:on nodes are adjusted, we explore how radius uncertainty propagates to fluid dynamics.To do so, we fit a normal distribu:on (,  %̂" ) to the op:mal region represen:ng the vessel radius (Figure 6).Sampling from a normal distribu:on is jus:fied by systema:cally analyzing PDF in the op:mal region.In vessels with reassigned radii, we sample from the distribu:on (,  %" ) , where  = ̂/̂$ is a parentdaughter scaling factor.Within the normal distribu:on that represents possible radii values, we sample 1000 :mes in the chosen representa:ve network and 100 :mes in other segmenta:ons.

RESULTS
Varia9on in 3D Rendering.One-way ANOVA (Table 2) of the radius and length obtained from the raw networks generated from VMTK (before junc:on and radius adjustment) show that the variability between the dimensions generated from each segmenta:on is more significant than what would be expected by solely random chance (Kauffman and Schering, 2007).The ANOVA analysis on seven pulmonary arteries matched in each segmenta:on shows sta:s:cally significant differences in length and radii.In the aor:c network, the radius of the innominate, brachial, right common caro:d (RCC), and right subclavian arteries differ significantly.The length of the leU common caro:d (LCC), brachial, leU vertebral, and RCC also varies significantly between segmenta:ons.This analysis demonstrates that the raw data obtained from VMTK is uncertain despite having the same user segment all images using standard segmenta:on parameters, methods, and soUware.This jus:fies the need to generate a method to extract length and radii values for vessels in a network while accoun:ng for uncertainty.Main (MPA), right (RPA) and leU (LPA) pulmonary arteries; Right (RIA) and leU (LIA) interlobular arteries; and right (RTA) and leU (LTA) trunk arteries.
Network Pruning. Figure 7 shows the impact of standardizing the network size between segmenta:ons.Results show predic:ons of pressure and flow at the midpoint of the MPA, RPA, and LPA.Before pruning, predicted pressures vary significantly between segmenta:ons with systolic differences of over 15 mmHg (Figure 7, solid lines).The segmenta:on with the highest number of vessels pre-pruning has a pressure range between 18 − 31 mmHg.The segmenta:on with the smallest number of vessels has a pressure range of 5 − 15 mmHg.Once all vessels are pruned to 167 vessels, the pressure range between each segmenta:on lessened to about 5 mmHg (Figure 7, dashed lines).Since a fixed flow profile is imposed at the inlet to the network, there is significantly less varia:on in flow than pressure.However, pruning s:ll impacts flow predic:ons favorably.Before pruning, the peak RPA flow varies by approximately 100 mL/s and aUer pruning, the varia:on becomes less than 50 mL/s.Flow discrepancies in the MPA and LPA flow are smaller.These results are summarized in Table 3.
Variance in VMTK Radius Es9mates.Figure 8 shows the CV and KDE for the pulmonary and aor:c vasculatures.The largest vessels in the pulmonary vasculature, with radii between 1−1.5 cm, have rela:vely small CVs, indica:ng that the centerlines generated from the segmenta:on and radius extrac:on for these vessels have less noise.For vessels with radii less than 0.5 cm, CVs range from 0−0.6, with one outlier reaching a CV of over 1.This vessel has a standard devia:on greater than its mean, which can occur when terminal vessels are close to image resolu:on limits.Smaller vessels consistently have a higher and broader range of CV values in each segmenta:on.
In the aor:c vasculature, CVs range from 0 − 0.32.Vessels with radii less than 0.4 cm and greater than 1 cm have a CV of less than 0.18.Vessels with radii between 0.4−1 cm have CV values greater than 0.27.Similar to the pulmonary arteries, the largest aor:c vessels have the smallest CV of less than 0.1.The pulmonary KDE plots have a unimodal, right-skewed distribu:on with the majority of vessel radii from all segmenta:ons being less than 0.5 cm, peaking at 0.1 cm.Given the varia:ons in segmenta:ons, the radii distribu:on varies across networks even aUer pruning.For example, in one segmenta:on, a radius of 0.1 cm is observed in 17.5% of vessels, while in another, this radii value is found in 12.5%.However, it should be noted that prior to pruning this difference more pronounced with up to 20% varia:on between networks (Figure 3).However, it should be noted that prior to pruning this difference more pronounced with up to 20% varia:on between networks (Figure 3).For the aor:c networks, the KDE follows a right-skewed bimodal distribu:on with a major and minor modes at 0.25 and 0.6 cm, respec:vely.Most aor:c vessels have a radius of less than 0.5 cm, with another large por:on having a radius between 0.5 − 0.75 cm.Similar to the pulmonary vessels, dis:nct networks have different radii distribu:ons.Despite that, the radii distribu:on between segmenta:ons is less than 5%.Junc9on Node Varia9on.Pressure and flow waveforms calculated at vessel midpoints before (solid lines) and aUer (dashed lines) modifying junc:ons are shown in Figure 9.The varia:on before and aUer modifying junc:ons is significantly smaller in the pulmonary vessels compared to the aorta.The pressure varies by 1 − 3 mmHg in the MPA, RPA, and LPA.The difference between flow predic:ons is also negligible in the pulmonary compared to the aor:c vasculature, but there is a dis:nct effect on the shape of the flow predic:ons.Before modifying junc:ons, pulse pressure in the aor:c vasculature is about 80 mmHg, ranging from 40 − 120 mmHg.AUer this change, the pulse pressure is 50 mmHg, ranging from 50 − 110 mmHg.Before adjus:ng the junc:ons, the flow to the LCC and brachial arteries are higher, whereas, aUer this modifica:on, more flow is directed to the descending aorta.Table 4 lists the range, minimum, and maximum pressure, and flow, before and aUer adjus:ng the junc:ons.Wave reflec:ons are also impacted by modifying junc:on nodes, par:cularly in the LCC.This is likely a result of reduced flow to the LCC aUer the junc:on modifica:on.This difference is amplified in the aor:c segmenta:on reported with a solid black line with a sizeable nega:ve dip and one large reflec:on around 0.6 s.AUer modifying junc:ons, we observe one steady reflec:on between 0.5 − 0.7 s.The opposite trend is apparent in the brachial artery, where predic:ons before junc:on modifica:on show one reflec:on and those aUer show two minor reflec:ons.Change Point Analysis.Figures 5 and 6 depict the results of the change point analysis used to iden:fy vessel radii.We choose a maximum of three change points in each vessel to balance model complexity and interpretability.We analyze nearly 1000 vessels within five segmenta:ons of pulmonary arteries (167 vessels each) and ten of the aortas (13 vessels each).We find that 9% of pulmonary arteries and 33% of aor:c vessels have either zero, one, or two change points.The rest of the vessels have three change points.AUer combining change point detec:on with our radius extrac:on protocol, less than 1% of vessels require post-processing during the network consistency check.This demonstrates the robust nature of the framework we have developed.Analysis of Raw Networks from VMTK.Analysis of the raw networks generated by VMTK from the 3D rendered surfaces shows that the CV is higher in the pulmonary vasculature compared to the aor:c, especially in the smallest vessels (Figure 8).This is an expected finding since the aor:c network only consists of rela:vely large vessels, while the pulmonary network has many small vessels near image resolu:on.Moreover, high devia:on is expected since pulmonary tree segmenta:on requires manual components.
The bimodal nature of the KDE plots in aor:c vessels (Figure 8) stems from the change in the size of the vessels analyzed.The ascending aorta, arch I, arch II, and descending aorta have significantly larger radii than the branching head and neck vessels.In contrast, the pulmonary vasculature has a rapidly branching structure, giving rise to a unimodal right-skewed KDE.These findings reflect uncertainty in the radii extracted from centerlines, obtained from the segmenta:ons of MRI and CT images.The importance of quan:fying uncertainty associated with quan::es extracted from image segmenta:on through the use of KDEs and PDFs has been supported by several studies.Baumgartner et al. (2019) generate PDFs for each segmenta:on and compare it to the distribu:on of the actual image.They argue that knowing these PDFs can greatly improve image segmenta:on accuracy (Baumgartner et al., 2019).KDEs have also been used for machine-learning applica:ons of image segmenta:on.Wang et al. (2013) leverage KDEs of abdominal blood vessels from CT angiograms to learn target distribu:ons and es:mate vessel loca:ons in segmenta-:ons of the same vessels in different pa:ents.The present study uses KDEs to understand vessel distribu-:ons within and between segmenta:ons, and we note that this analysis can be a powerful tool for gaining informa:on about vessel uncertainty.
It is well known that limited image resolu:on, segmenta:on inaccuracies, or inexactness in centerline genera:on (Lins et al., 2020;Sharma and Aggarwal, 2010;van Rikxoort and van Ginneken, 2013;Buelow et al., 2005) impact 3D rendering.Other factors, including noise, poor contrast, and blurred boundaries, also increase the difficulty of extrac:ng the smallest vessels (Cui et al., 2019;Lesage et al., 2009;Hirsch et al., 2012).In the study by Colebank et al. (2019) examining image segmenta:on uncertainty, the uncertanty increases as vessel radius decreases, sugges:ng that small vessels are more sensi:ve to segmenta:on parameters.Van Horssen et al. (2016) postulate that small vessels are more likely underes:mated due to image resolu:on.Schwarz et al. (2020) found that in the human coronary arterial circula:on, the smallest vessels are subject to segmenta:on errors, resul:ng in triangular loops and spurious side branches that must be corrected in post-processing.These studies validate our findings that the aor:c networks with larger vessels have less varia:on than the small vessels in the rapidly branching pulmonary network.This varia:on is not a flaw of the soUware used to generate segmenta:ons and networks, but points to the need for post-processing when using image-generated networks for hemodynamic predic:ons.Network Standardiza9on.Standardizing networks to have a consistent number of vessels is impera:ve to enable comparison.AUer pruning pulmonary networks, we find that fluid dynamic differences stemming from segmenta:on varia:on become less apparent.An alterna:ve way to account for large vessel network size differences is to adjust structured tree parameters, α and β, genera:ng networks with similar total vessel count and a similar total vessel area or volume.However, this study inves:gates changes in geometric quan::es, such as radius and length, keeping parameters constant to emphasize the uncertainty caused by vascular geometry.Thus, we prune networks rather than adjus:ng α and β parameter values, though this can be inves:gated in future studies.Miller et al. (2023) argue that networks with disparate branch numbers are an imaging consequence rather than a biological factor, and thus, pruning is essen:al when comparing different segmenta:ons.In addi:on, they find that more informa:on about networks can be gained by applying pruning techniques (Miller et al., 2023).Colebank et al. (2019) find that choices in segmenta:on parameters greatly impact the number of vessels in the network and that network size affects the model output.Since mean pressure and flow metrics are used to diagnose diseases, such as pulmonary hypertension (Hoeper et al., 2013), standardizing results and assessing uncertainty is cri:cal.This also highlights the importance of determining the minimum number of imaged arteries needed for reliable diagnos:c results.If a subset of larger vessels can yield valid, comparable results, excluding smaller vessels that exhibit dras:c imaged-related uncertain:es may be beneficial.This idea has been studied by several groups, who use surrogate or reduced order models to lessen computa:onal complexity (Paun et al., 2021;Borowska et al., 2023;Ragkousis et al., 2016;Epstein et al., 2015).
Table 5: Range, minimums, and maximums of pressure and flow predicXons comparing the "true" geometry to pre-dicXons that are from sampled geometries.The "true" geometry comes from the mean of the opXmal secXon selected by change points and the sampled comes from sampling from a normal distribuXon around the mean.
Junc9on Node Placement.The placement of nodes within junc:ons impacts both the length and radii of vessels.Large vessels, especially the aorta, are known to taper along their length (Caro et al., 1978).Due to the rapid change in vessel radii, junc:on node placement in these vessels is essen:al.Our algorithm addresses this by shiUing the VMTK-derived nodes to the center of the os:um.As a result, sec:ons of centerline data previously characterized as part of daughter vessels become reclassified as part of the parent.Modifying junc:on nodes impact hemodynamics in the pulmonary and aor:c vasculature (Figure 8) but impacts aor:c vessels more significantly.Since aor:c vessels are longer and wider than pulmonary vessels, geometric changes resul:ng from modifying junc:ons are more pronounced, and this effect is propagated to hemodynamics.In par:cular, in the descending aorta, one segmenta:on has a nearly 60 mmHg varia:on in systole between junc:on node loca:ons.This is not the case in pulmonary vessels, with an average pressure varia:on of only 3 mmHg in systole.The wave reflec:ons in the aor:c flow predic:ons are also affected.Abdullateef et al. (2020)  reflec:ons, finding that as vessel radius decreases, wave reflec:ons increase and begin to overlap.Accurately placing nodes in tapering vessels is crucial, as an incorrect placement can cause wrong reflec:ve wave proper:es.Due to the well-known interplay between geometry and hemodynamics, genera:ng accurate vessel geometry is impera:ve.Although this study uses the cylindrical Navier−Stokes equa:ons to model hemodynamics, considering the concept of steady flow within a cylinder, namely Poiseuille's law q = πpr 4 /8µl (Pfitzner, 1976;Demers and Wachs, 2022), allows us to gain insight into how radius and length influence pressure and flow dynamics.Clearly, q and p are directly impacted by r and l in Poiseuille's law.Notably, r is raised to the fourth power, whereas l is raised to the first power.These results are reflected in our study, though it should be noted that changing the junc:on nodes' loca:on affects the vessel length and radius since nodes with smaller or larger radii may be moved to a different part of the vasculature.The laMer, in par:cular, affects the aor:c vasculature, where small and large vessels merge at junc:ons.
Change Point Analysis.We sequen:ally combine two methods for change point analysis, u:lizing R/changepoint to detect the op:mal number of change points (Killick and Eckley, 2014) and R/segmented to place them (Muggeo, 2008).Both methods are computa:onally fast, and their combina:on enables accurate and automated assessment of vascular data.R/changepoint determines the op:mal change point count using mean and variance to detect changes in the data.However, it is an approxima:on algorithm, only giving an es:mated change point placement rather than an exact loca:on (Killick and Eckley, 2014).In contrast, R/segmented accurately places change points through a customizable and flexible algorithm, where users define linear regression models based on prior knowledge of data structure (Muggeo, 2008).This package requires users to input the number of change points detected within each dataset (Muggeo, 2008), so we use the result from R/changepoint as an input.Coupling these two packages allows us to create an automated pipeline for accurate vascular network analysis requiring minimal user input.This is advantageous when comparing mul:ple vascular datasets, which may encompass thousands of vessels.
By limi:ng the number of change points to three, we ensure that detected points capture meaningful transi:ons in a vessel's structure, such as the os:um region or image resolu:on loss, while preven:ng the overfikng of the data.It enables consistent and computa:onally fast analysis, allowing for the development of a robust algorithm for extrac:ng radii.While we could define exactly three change points in all datasets, this leads to inaccurate change point placement in vessels with an op:mal count of less than three, as observed in 9% of pulmonary vessels and 33% of aor:c vessels.
Radius Sampling.Figure 10 shows that vessel radius greatly impacts hemodynamics, following the paMern we expect from Poiseuille's law.Given the predicted pressure and flow values for various radii, a small change in radius causes a large change in hemodynamics.For example, when we perturb the LCC's radius following a distribu:on N(0.36,0.01),we observe a 70 mmHg range in systolic pressure.This is also shown in a study by Colebank et al. (2019), where even slight changes in vessel dimensions influence pressure and flow.Secomb (2016) notes that small changes in diameter can widely modulate blood flow and emphasizes the need to precisely iden:fy vessel dimensions to understand how blood is distributed throughout the body.We note that pulmonary blood pressures fluctuate throughout the respiratory cycle, evident when using right-heart catheteriza:on for pulmonary artery measurements (Kovacs et al., 2014).Pressure measurements are greatly impacted by pa:ent angle during data collec:on, changes in al:tude, ambula-:on, and pre-exis:ng heart arrhythmias.Especially when measuring the same subject over :me, this can lead to a pressure varia:on of between 5−15 mmHg (Magder, 2018;Alam et al., 2022).Our pressure range obtained from sampling radii falls within the an:cipated range, depending on how much the radius is changed.It is common prac:ce to segment cardiovascular networks to build pa:ent-specific models used as noninvasive methods to monitor and understand disease (Colebank et al., 2021a;Taylor-LaPole et al., 2023;An:ga et al., 2008;Carson et al., 2019;Dobroserdova et al., 2016).Our findings demonstrate that accurate radius measurements are necessary to tune models to pa:ent data accurately, especially when using pa:ent-specific parameters, such as vessel s:ffness.
The use of change points to iden:fy the vessel radius and its standard devia:on is novel, providing an automated way to iden:fy the vessel radius accurately.Most previous studies either average all nodes within the vessel (Pfaller et al., 2022;Mulder et al., 2011), or assign the vessel radius using a set propor:on of the nodes within the vessel without accoun:ng for different vessel structures (Colebank et al., 2021b).
Correct iden:fica:on of vessel radii is important beyond computa:onal modeling.In the current study, our algorithm is used for iden:fying vessel radii in vasculature networks informing the geometric domain for computa:onal analysis.However, many studies segment vascular structures from CT and MRI images without using CFD.This is useful when studying diseases and surgical interven:ons, such as hypoplas:c leU heart syndrome (HLHS) and chronic thromboembolic pulmonary hypertension (CTEPH).HLHS pa:ents oUen have enlarged and/or deformed aortas due to reconstruc:ve surgery (Taylor-LaPole et al., 2023).Our algorithm can iden:fy pa:ents with abnormal vessel radii and compare radius distribu:on to control individuals, such as DORV pa:ents.Another example is quan:fying vascular remodeling in pa:ents with CTEPH, including quan:fying enlargement of the MPA an acute radius reduc:on in vessels with ring-or web-like lesions.Compu:ng radii values with uncertainty can provide valuable guidance for surgical interven:ons, assis:ng physicians in iden:fying which vessels with lesions should be targeted.
Limita9ons and Future Studies.This paper outlines a semi-automa:c framework to extract and refine pa:ent-specific geometries from medical images.We use a linear regression model to define change points for all vessels.With these change points, we can avoid unclear os:um regions and reliably detect representa:ve radii.However, since we use an exponen:al func:on to model tapering radii, we can reduce computa:onal complexity by directly defining change points with an exponen:al model in those cases (Chen and Gupta, 2011).
In addi:on, it must be noted that we make several assump:ons regarding the evolu:on of a vessel's radius and the loca:ons of junc:ons to develop this framework.These assump:ons are made based on a thorough analysis of typical vessel structures in large networks and a literature review.However, our framework is highly customizable, so these assump:ons can be adjusted to suit the needs of various vascular networks that differ from the aorta and pulmonary arteries, such as the renal, cerebral, and hepa:c arteries.Future studies will inves:gate the impact of specific modeling assump:ons on hemodynamic predic-:ons in different systems.
Another limita:on is that we sample radii for each vessel independently.While this choice allows us to inves:gate how geometry impacts hemodynamics, it overlooks the in:mate connec:on between vessels in the system.When we sample independently, daughter radii can exceed their parent's, viola:ng the morphometric condi:on.In the future, we will account for covariance enabling us to consider the link between parents and daughter vessels when sampling (McNeil, 2008).
We generate a rela:vely small number of segmenta:ons for analysis in this study due to the :me-consuming nature of manual segmenta:on of complex structures.We plan to create more segmenta:ons of pulmonary and aor:c vessels to test our method further.To do so, we will use automa:c segmenta:on soUware, such as (Polek et al., 2022), and compare this to the manual segmenta:ons we have already completed.Another component is studying how changes in segmenta:on threshold and smoothing impact predic:ons.This was analyzed by Colebank et al. (2019) in excised mouse pulmonary artery data captured by microcomputed tomography images.Colebank's results had similar varia:on in flow and pressure compared to those reported in this study.Finally, it is important to factor in the noise in images.Animal dta are oUen obtained under anesthesia, while human data are obtained in the awake state.The laMer may have higher impact on image resolu:on than analysis of radii extracted from the images.

CONCLUSION
This study develops a robust framework to generate pa:ent-specific labeled directed trees from medical images and obtain vessel measurements for use in fluid dynamics simula:ons while considering uncertainty.We find notable varia:on in the raw networks generated using VMTK, especially in the smallest vessels, when segmen:ng the same images mul:ple :mes despite using standard segmenta:on procedures and having the same user segment the images.To minimize this uncertainty, we develop a robust pipeline to prune networks, shiU junc:on nodes to the center of the os:um, and determine vessel radius using change points.Our results emphasize the significance of post-processing and quan:fying uncertainty when genera:ng networks for medical applica:ons.Current diagnos:c procedures for many diseases require invasive protocol, including right-heart catheteriza:on.This study demonstrates that careful analysis of image analysis data combined with CFD modeling has the poten:al to augment or mi:gate the need for invasive studies via in-silico simula:ons.Despite the poten:al clinical advances CFD provides, medical conclusions should only be made by first assessing the level of uncertainty present in diagnos:c informa:on, fully understanding limita:ons and possible errors.In the absence of fixed geometric values of vessels due to scanning ambiguity, researchers must navigate uncertain measurements by developing effec:ve methodologies, such as the one we have outlined here, for increasingly accurate analyses.

Figure 1 :
Figure 1: Workflow for our methodology, consisting of image segmentation and analysis, creation of labeled directed tree from centerlines, junction node modification, change point analysis, radius identification, and construction of a computational domain for fluid dynamics simulations.

Figure 2 :
Figure 2: 3D rendering from (a) a normotensive pulmonary arterial network and (b) the aortic vasculature of a double outlet right ventricle patient.Centerlines (black lines) and structured trees, used outside the imaged region, constitute a labeled directed tree.

Figure 3 :
Figure 3: Kernel density estimates in five pulmonary arterial networks for trees (a) before and (b) after pruning to have the same number of vessels.
2: Star:ng at the first node along each daughter vessel  ,! , we iden:fy all nodes  = 1, … ,  for which  *  &/0 ⁄ ≤ 0.1.Here,  denotes the index of first the node for which  *  &/0 ⁄ > 0.1.The new nodes along the parent branch are computed by averaging the loca:on of the daughter vessel coordinates

Figure 4 :
Figure 4: Junction node adjustment in 3D space for a bifurcating vessel.Dark gray lines and teal circle represent the VMTK-generated centerlines and junction node  !, respectively.(a) Nodes along daughter vessels, where daughter vessel 1,  " , is in red and daughter vessel 2,  # , is in blue, for which  $  %&' ⁄ ≤ 0.1,  = 1, … , .(b) The green line and purple circle denote the new parent vessel and junction node  .! , calculated by an average between daughter

1 :
Sta9s9cal change point calcula9on.A change point  ∈ {1, … ,  − 1} can be iden:fied within the ordered sequence of radii data (determined from VMTK values)  = { !, … ,  4 } for which sta:s:cal proper:es \ !, … ,  5 ] and \ 56! , … ,  4 ] differ significantly.The dataset may contain mul:ple change points, represented as  = { !,  " , … ,  7 }, where  is the number of change points.We use a two-step process to determine the change points using publicly available R packages.We use R/changepoint to iden:fy the number of change points and R/segmented to define the loca:on of change points within each vessel.The two-step process is needed since R/segmented requires user-specifica:on of the number of change points, and R/changepoint does not precisely indicate where change points occur.

Figure 5 :
Figure 5: Binary segmentation is used to place change points along radii data attained from VMTK (grey points).In each panel, the blue lines are the piecewise linear fit determined in Equation (9), and the red circles denote the change point placement.Panel 1 shows the radii before change points are placed, and panel 2 shows the segmentation after placing the first change point.Panel 3 shows the placement of the second change point.It is placed to the left of the first change point as this placement leads to a higher maximum likelihood.Finally, panel 4 shows the optimal placement of three change points.

Figure 6 :
Figure 6: Radius extracXon for (a) pulmonary and (b) aorXc arteries.Change points (red circles) denote locaXons where the vessel radius changes significantly.The segment used to determine the vessel radius (green) avoids the juncXon region (pink box) and rapidly changing radii at the end of the vessel.In the aorta, esXmated inlet and outlet

Figure 7 :
Figure7: Pressure and flow dynamics over a cardiac cycle for each segmentation at the midpoint of the MPA, RPA, and LPA in the pulmonary arteries before (solid lines) and after (dashed lines) standardizing network size through pruning.Each color represents a distinct segmentation.

Figure 8 :
Figure 8: Coefficient of variation and kernel density estimates for (a) pulmonary arteries and (b) aortic vessels.Each color represents a distinct segmentation.

Figure 9 :
Figure9: Pressure and flow dynamics over a cardiac cycle for each segmentation at the midpoint of (a) the MPA, RPA, and LPA in the pulmonary arteries and (b) the descending aorta, LCC, and brachial artery before (solid lines) and after (dashed lines) modifying VMTK-derived junctions.Each color represents a distinct segmentation.

Figure 10 :
Figure 10: Pressure and flow dynamics over a cardiac cycle at the midpoint of the (a) MPA, RPA, and LPA in the pulmonary arteries and (b) descending aorta, LCC, and brachial artery.Gray lines (1000) denote simulaXons with radii sampled from a normal distribuXon based on the mean and standard deviaXon of each vessel's extracted radius, and the colored line (1) represents the simulaXons with radii set at the mean radius based on our extracXon procedure.These simulaXons were based on a representaXve segmentaXon of the pulmonary and aorXc vasculatures.

Table 1 :
Parameters used in the fluid dynamics model for the pulmonary and aorXc networks.

Table 2 :
One-way ANOVA test results predict the error sum of squares (SSE) and p-values for seven vessels in the pulmonary vasculature and eight vessels from the aorXc network that can be matched in each segmentaXon.For conXnuity, all aorXc vessels (ascending, arch I, arch II, and descending) were treated as one exponenXally tapering vessel.P-values in bold text are staXsXcally significant under a 0.05 threshold.

Table 3 :
Range, minimum, and maximum pressure, and flow before and afer pruning.